Initialisation
## [1] "number of cores available = 10"
#Phi[1] ; eta = valeur de fin
#Phi[2] = valeur du noeud
#Phi[3] = echelle
m <- function(t, eta, phi) (phi[1] + eta)/(1+exp((phi[2]-t)/phi[3]))
#=======================================#
param <- list(sigma2 = 0.05,
rho2 = 0.1,
mu = c(5,4,0.5),
omega2 = c(0.5,0.1,0.01))
F. <- length(param$omega2) #dimension de phi
#=======================================#
t <- seq(2,6, length.out = 10) #value of times
# --- longitudinal data --- #
G = 40 ; ng = 4 #nombre de groupe et d'individu par groupe
n <- G*ng*length(t) ; N <- G*ng #nombre total de data, et #nombre total d'individu
dt_NLME <- NLME_obs(G, ng, t, param, m)
Y <- dt_NLME$obs

SAEM avec simulation par MCMC
\[\log f(Y, Z ; \theta) = \langle \Phi(\theta) ; S(Y,Z) \rangle - \psi(\theta)\]
#Petite fonction pour retourner rapidement l'appel Phi[attr(Phi, i)] où i est '1', '2', ...,
`%a%` <- function(x,var){
if(length(var)== 1) return(x[ attr(x,as.character(var)) ])
lapply(var, function(v) x%a%v) %>% unlist
}
Phi <- function(sigma2, rho2, mu, omega2, nu2)
{
p <-c(- n/(2*sigma2), -N/(2*rho2),
- G/(2*omega2), G*mu/omega2 )
return(p)
}
S <- function(eta, phi)
{
s <- c(mean((Y - get_obs(m, dt_NLME, eta = eta, phi = phi) )^2 ), #1
mean(eta^2), #2
apply(phi^2, 2, mean), apply(phi , 2, mean) ) #3 et 4
attr(s, '1') <- 1
attr(s, '2') <- 2
attr(s, '3') <- 2 + 1:ncol(phi)
attr(s, '4') <- 2 + ncol(phi) + 1:ncol(phi)
return(s)
}
Metropolis Hastings
loglik.phi <- function(x, eta, Phi)
{
s <- S(eta,x)
sum( (Phi*s)%a%c(1,3,4) ) #indice de S_1 et S_{3,.} puis S_{4,.}
}
loglik.eta <- function(x, phi, Phi)
{
s <- S(x, phi)
sum( (Phi*s)%a%c(1,2) )
}
SAEM
Initialisation
M <- 1 #nombre de simulation
u <- function(k) ifelse(k<200, 1, 1/(k-199) )
# --- Initialisation des paramètres --- #
para <- param %>% sapply(function(x) x + runif(1))
para$rho2 = 0.2 ; para$mu <- c(6,3,1) ; para$omega2 <- rep(.1,3)
# --- Initialisation des chaines MC : Z_0 ---
Z <- list(eta = 1:M %>% lapply(function(i) rnorm(G*ng, 0, param$rho2) %>% matrix(ncol = 1) ),
phi = 1:M %>% lapply(function(i) matrix(rnorm(F.*G, param$mu, param$omega2), nrow = F.) %>% t ) )
Boucle
sim <- function(niter, Phih, eta, phi)
{
M <- length(phi)
eta <- 1:M %>% lapply( function(i)
MH_High_Dim_para_future(niter, eta[[i]], sd = .036, loglik.eta, phi[[i]], Phih, cores = ncores-1 ))
phi <- 1:M %>% lapply( function(i)
MH_High_Dim_para_future(niter, phi[[i]], sd = c(.028, .036, .021), loglik.phi, eta[[i]], Phih, cores = ncores-1 ))
list(eta = eta, phi = phi)
}
maxi <- function(S)
{
list(sigma2 = S%a%1,
rho2 = S%a%2,
mu = S%a%4,
omega2 = S%a%3 - (S%a%4)^2 )
}
|
|
sigma2
|
rho2
|
mu1
|
mu2
|
mu3
|
omega21
|
omega22
|
omega23
|
|
Oracle
|
0.0510249
|
0.094816
|
4.97881
|
3.950752
|
0.4798243
|
0.5195325
|
0.0975481
|
0.0072938
|
|
Initialisation
|
0.5404334
|
0.200000
|
6.00000
|
3.000000
|
1.0000000
|
0.1000000
|
0.1000000
|
0.1000000
|
niter <- 50
MH.iter <- 10
res <- SAEM(niter, MH.iter, u, para, Phi, S, Z, sim, maxi, eps = 1e-3, verbatim = 2)
saveRDS(res, 'saem.rds')
gg <- plot(res)
## [1] "SAEM execution time = 15min 39sec"
Résultat de l’algo SAEM-MCMC
|
|
sigma2
|
rho2
|
mu.1
|
mu.2
|
mu.3
|
omega2.1
|
omega2.2
|
omega2.3
|
|
Valeur réelle
|
0.0500
|
0.1000
|
5.0000
|
4.0000
|
0.5000
|
0.5000
|
0.1000
|
0.0100
|
|
Valeur estimée
|
0.0525
|
0.2099
|
4.9736
|
3.9507
|
0.4825
|
0.2117
|
0.0945
|
0.0072
|
|
Rrmse
|
0.0491
|
1.0995
|
0.0053
|
0.0123
|
0.0350
|
0.5766
|
0.0548
|
0.2778
|


